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ESTIMATION  OF  TIME-OF-ARRIVAL 
OF  UNDERWATER  ACOUSTIC  SIGNALS 
BY  SPLINE  FUNCTIONS 


INTRODUCTION 

An  important  problem  in  underwater  signal  processing  research  is  the  determination 
of  the  signal  onset  (or  time-of-arrival)  of  an  acoustic  signal  from  possib'y  noisy  discrete 
data.  Once  an  accurate  estimation  of  this  value  is  obtained,  the  transducer  and  panel 
transfer  function  identification  techniques  can  be  applied.  In  the  report  [1],  the  underwa¬ 
ter  acoustic  signal  is  represented  by  a  spline  curve  with  equally  spaced  knots  such  that 
the  initial  knot  that  lies  in  the  interior  of  the  time  interval  clearly  indicates  when  the 
curve  “takes  off’  to  the  right,  and  hence,  can  be  used  to  define  the  signal  onset.  In  par¬ 
ticular,  if  the  signal  arrives  fairly  sharply  a  linear  spline  curve  can  be  used  to  model  the 
signal.  However,  if  the  signal  arrives  and  increases  in  magnitude  smoothly  and  slowly,  a 
higher  order  spline  curve  such  as  a  cubic  spline  is  recommended.  In  a  recent  report  [2], 
it  was  pointed  out  with  several  interesting  examples  that,  indeed,  a  cubic  spline  curve 
does  not  give  an  accurate  estimation  of  the  signal  onset  when  the  signal  function  rep¬ 
resentation  has  nonzero  first  or  second  derivatives  at  the  time-of-  arrival.  Sir  e  a  cu¬ 
bic  spline  curve  provides  a  better  model  than  its  linear  counterpart  to  curve  fitting,  and 
hence,  to  representing  an  acoustic  signal,  we  will  modify  the  cubic  spline  model  in  Ref. 

1  by  using  “stacked  knots”  at  the  initial  knot  position.  As  a  consequence,  the  B-spline 
series  is  increased  by  two  terms  which  may  both  drop  out  when  the  arriving  signal  has 
zero  first  and  second  derivatives  at  the  signal  onset,  so  that  the  approximation  procedure 
is  adaptive  in  nature. 

The  main  objective  of  this  report  is  to  provide  useful  algorithms  for  determining 
the  signal  onset  from  discrete  measurement  of  the  acoustic  signal  which  is  contaminated 
with  noise.  The  method  of  penalized  least-squares  discussed  in  the  report  [1]  will  be 
studied  in  great  details.  In  particular,  all  relevant  matrices  for  the  setting  of  this  partic¬ 
ular  problem  will  be  analyzed  so  that  the  method  of  generalized  cross-validation  (GCV), 
(See  [3-5]),  can  be  modified  and  extended  for  our  study.  A  global  search  procedure  will 
then  be  used  to  determine  a  good  estimation  of  the  signal  onset.  For  the  sake  of  com¬ 
parisons,  algorithms  based  on  SVD  (singular  valued  decomposition)  and  tridiagonization 
will  both  be  proposed.  In  addition,  whenever  it  is  appropriate,  flow  charts  will  be  given 
to  illustrate  the  algorithms. 

SPLINE  REPRESENTATIONS  OF  ACOUSTIC  SIGNALS 

Let  [0,d]  denote  the  time  interval  on  which  an  underwater  acoustic  signal  is  mea¬ 
sured,  and  let  c  =  fo-0  <  c  <  d,  represent  the  signal  onset  tor  time-of-arrival)  which  is 
to  be  determined.  Then  the  acoustic  signal  is  represented  by  a  spline  curve  with  initial 
knot  at  to-  The  knot  sequence  of  this  spline  curve  was  given  on  page  7  of  Ref.  1  by 

th'-  to  <  *i  <  •  •  •  <  tn  =  d  <  t„  + 1  <  •  <  t„  +  k 

where  0  <  f0  =  e  <  d.  Note  that  k  knots  tn+I .  •  •  •  ,  tri  +  i  are  tacked  on  to  the  right 
of  the  measurement  interval  [0,rf]  to  allow  us  to  represent  the  last  k  B-splines  in  the 


1 


5-spline  series  representation  of  the  spline  curve  with  degree  k.  Recall  that  a  5-spline 
with  degree  k  is  supported  on  k  +  2  knots,  so  that  the  last  5-spline  has  support  given 
by  n+t]  which  overlaps  with  the  time  interval  [0, d\.  As  suggested  by  George  and 

Muise  [2],  it  is  quite  possible  that  the  arrival  of  an  underwater  acoustic  signal  is  not  at 

all  smooth  in  the  sense  that  some  of  the  initial  derivative  values,  /'(to)- 1  •  •  o). 

of  the  function  /(/)  which  represents  the  signal  do  not  vanish.  Hence,  in  order  to  use  a 
A’-th  degree  spline  curve  to  represent  the  acoustic  signal,  we  recommend  in  this  report  to 
extend  the  5-spline  series  to  the  left  by  stacking  k  —  1  knots  at  the  initial  knot  so  that, 
on  one  hand,  non-zero  initial  derivative  values  /o)(/0),  1  <  j  <  k  —  1,  are  allowed,  and 
on  the  other  hand,  the  spline  curve  still  takes  off  at  the  initial  knot  t0.  That  is.  we  now 
consider  the  knot  sequence 

th:  t-k+i  =  •••  =  <-,  =t0  <  ti  <•••<*„  <  •  <  /„+*.-  (1) 

where,  again,  0  <  tQ  =  c  <  d  and  tn  =  d.  For  the  distinct  knots  f0.  ■  •  •  .t„+k,  we  will  still 
assume  that  they  arc  equally  spaced  as  in  Ref.  1;  that  is,  for  j  =  1.  •••.?;  +  k,  we  set 

,  d  —  c 

tj  —  tj  —  \  +  h,  h  —  - , 

n 

since  there  is  nothing  to  gain  by  making  the  problem  more  complicated. 

B-Splines  With  Non-Zero  Initial  Derivatives 

To  construct  the  5-splines  with  the  new  knot  sequence  t /, ,  we  return  to  Eq.  (9)  on 
page  5  of  Ref.  1  by  setting 

Bk.j(t)  =  Bk'thj(t)  =  Nkfyt  ~  c)  (2) 

where  j  =  0,  •  •  • ,  n  —  1  and  Nk ( t )  denotes  the  5-spline  with  degree  k  and  knots  at  the 
integers  having  support  given  by  the  interval  [0,  k  +  1],  (See  Fig.  1  below  and  on  page  5 
of  Ref.  1 ). 


Fig.  1  -  5-Splines  of  orders  1  -  4 

For  the  knots  which  are  no  longer  equally  spaced,  we  must  compute 

the  5  splines  Bk,^(t).  J  =  —k  4-  without  relying  on  A\fM.  For  this  purpose, 

the  method  described  on  page  12  of  Ref.  G  works  quite  well,  since  Bernstein  coefficients 
must  be  obtained  for  later  purposes.  In  this  report,  we  will  only  consider  linear  and  cu¬ 
bic  splines  (i.e..  k  —  1  and  3.  respectively).  Since  no  stacked  knot'-  are  necessary  for 
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linear  splines,  only  cubic  B-splines  B3j(t),j  =  -2,-1,  will  be  given.  [For  convenience, 
we  will  use  the  notation  B}(t)  =  B3 .>(/)-] 

Using  Bernstein  representations  as  in  Ref.  1,  we  have  (for  k  =  3): 

B-2(t)  =  =  Jv3Q(t  -  c>)  (3) 

and 

B-i(t)  =  B3.t*,_i(0  =  -  c)),  (4) 

where  the  Bernstein  coefficients  of  Ar3  and  N3  are  given  in  Fig.  2. 

N3  N  3 


0  1  \  A  0  0  0 

1 - 1 - —1 


0  0  i 
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I  I  I  0  0  0 

_ I _ _ I 


0  1  2 


0  1 


2  3 


Fig.  2  -  Bernstein  coefficients  of  N3  and  N3 
Of  course,  as  mentioned  in  Eq.  (2),  for  j  —  0, 1,  •  •  • ,  n  —  1,  we  have 

Bj( t)  =  B, ,.»,>(«)  =  JV3  (i(t  -  c)  -  j),  (5) 

and  the  Bernstein  representation  of  N3(t)  was  already  computed  on  page  6  of  Ref.  1  as 
shown  in  Fig.  3  below. 


0  0  0  A  | 

i _ I _ 
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N3(x) 

Fig.  3  -  Bernstein  representation  of  A r3(t) 
In  the  usual  Cartesian  formulation,  we  may  also  write: 


-v.(0 


3(1  -  t)2t  +  |(1  -  t)t2  + 

4(2 -n3 

0 


if 

if 

otherwise 


0  <  t  <  1 
1  <  t  <  2 


3 


Chui 


()</<! 

1  <  t  <  2 

2  <  t  <  3 


and 


-V:i(0 


l  0  otherwise. 

Lsing  the  notation  in  Eqs.  (3)  to  (5),  any  cubic  spline  function  S^(t)  that  represents  the 
underwater  acoustic  signal  f(t)  with  signal  onset  at  t0  =  c  and  possibly  nonzero  values 
of  f'(t o)  and  /"(to)  is  given  by  the  cubic  spline  series 

71-1 

Ss(t)  =  cJDJ{t]  (C) 

j=-  2 

with  knot  sequence 

where  0  <  /„  =  c  <  <7.  tn  =  d,  and  h  =  (d  —  c)/n.  Observe  that  Eq.  (C)  differs  from 
Eq.  (11)  on  page  7  of  Ref.  1  in  that  two  extra  terms  have  been  introduced  to  allow  an 
adaptive  curve  fitting  scheme  by  this  new  cubic  spline  representation. 

The  Coefficient  Matrices  in  L 1  Approximation 

For  L “  =  L2  0.  d]  least-squares  estimation  of  underwater  acoustic  signals,  it  is  essen¬ 
tial  to  compute  the  coefficient  matrices 

AU'  =  [?q;] 

where  —  k  +  1  <  j  <  n  —  1  and 

tfj—  f  Bk.th  ,j(  1  v  t  •  :si 

,/u 

Hence.  for  linear  sp] i i ios ,  where  k  —  1.  the  coefficient  matrix 


0  otherwise 


i  f3 

6 

if 

0  <  /  <  1 

|(2  —  t)3  T(2  —  t)2(i  —  l)  +  2(2  —  t)(t  — 

l)-  +  l(t-l)3 

if 

1  <  f  <  2 

f  (3  —  /)3  +  2(3  -  t)2(t  -  2)  +  (3  -  t)(t  - 

2)2  + i(/-2)3 

if 

2  <  t  <  3 

if 

3  <  t  <  4 

|(i  +  nt3 


-^3  ( f ) 


17J  ( —  7 ) 3  +  2(2  —  f)2(/  —  1)  +  (2  —  /)(f  —  1)-  +  ~(t  —  1 ) 3  if 

i(3-/)3  if 
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■4  1  . 

1  4  .  '  • 


.4,,«  =  5 


.  4  1 

'  1  2. 


given  in  Eq.  (39)  on  page  15  in  Ref.  1  is  unchanged.  For  cubic  splines,  where  k  =  3. 
however,  the  dimension  of  the  matrix  /I3  ^  is  increased  by  2,  becoming  an  (/;  +2)  x  ( 11  +  2 ) 
matrix.  In  the  next  subsection,  we  will  snow  that 


[compare  with  Eq.  (40)  on  page  15  of  Ref.  l],  where 

p  _  [2232  1575 

1  *  «-r  r*  r»r\A  t 


is  a  2  x  2  block, 


is  a  2  x  ( n  —  3)  block, 


£b  = 


Pi 

0 

i»T 

Cn-3 

D 

0 

dt 

e3 

f.  l],  where 

'2232 

1575' 

1575 

3294 

3 

0  0 

... 

239 

2  0 

2416 

1191 

120  1 

1191 

2416 
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120 

1191 

•  .  * 
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•  •  • 
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•  • 

*  *  *  • 

120 

•  •  « 

1191 

0 

• 

1  120  1191 

2416 

is  an  (n  —  3)  x  (n  —  3)  banded  Toplitz  symmetric  matrix  as  Eq.  (41 )  on  page  15  of  Ref. 

1. 


D  =  1  0 

120  1  0 

.1191  120  1. 

an  in  -  3)  X  3  Toeplitz  matrix  as  in  Eq.  (42)  on  page  15  of  Ref.  1.  an 


2396 

1062 

60 

1062 

1208 

129 

60 

129 

20 

o 


C’hui 


is  a  3  X  3  block  as  in  Eq.  (43)  on  page  15  of  Ref.  1.  Of  course,  D[  denotes  the  transpose 
of  D\  in  Eq.  (12)  and  Dr  the  transpose  of  D  in  Eq.  (14). 

Derivation  of  the  Coefficient  Matrix 

To  derive  the  coefficient  matrix  in  Eq.  (10),  the  Bernstein  representations  of 
A\(O.AT3(<)  •  and  A’3(/)  shown  in  Figs.  2  and  3  are  used  together  with  the  relationships 
in  Eqs.  (3)  to  (5)  in  the  integral  in  Eq.  (8)  for  b3}  (with  k  =  3).  After  a  change  of  vari¬ 
ables  from  B3jh  l(t)  to  iV3 (  / ) ,  .Y3 (f  ),  or  .V3 (t),  the  integrals  can  be  evaluated  by  using  a 
formula  in  Eq.  (37)  on  page  12  in  Ref.  1,  namely: 


f  '  +  l  p  n  _  (^’')2  (l  +  p)'(m  +  ?)■  k  ,k 

Jr  kQk  (2A-  +  1  £\m\p\q\  tm  " 


£-\-m~k  p-\-q~k 


[£  +  p)\{m  +  g)\  k  k 
(\m\p\q\  C(m  pr 


where  k  =  3  is  used  for  cubic  splines.  In  the  following,  the  superscript  k  for  C(m  and 
will  be  dropped  for  convenience. 

(a)  For  b3_2  _2,  we  use  N3(t)  so  that 


yielding: 


k,m}  =  (o,i.i,l),  {<*„}  =  (i.o, 0.0). 


>--  =  ^[20.(l)%12  +  0.(i)+4(I) 
+  9,0  +  12.(^  +  lo. (>).(!) 

+  4.(i)  +  10.(l)  (I)+20.(i)2!/, 


=  ^  •  2232. 


(b)  For  ft'l2  _i  =  t>3_ i  _2 ,  we  use  both  1V3  ( f )  an<l  A 7j(f)  so  that  in  addition  to  Eq.  (17). 
we  also  need 


(eon)  -  (0.0. ).  {,/  }  =-  (— 

*  ’  V  2  12/  ’  V 12  G  G  G 


yielding: 


G 
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&3-2,-l  =  ^>-1,-2  = 


•{20(co3Co3  4  do 3 do 3 )  4-  10(ei2C03  -f  dj^do 


4  4(c2iCo3  4  d2i<io3)  4  (030^03  +  dsodoz)  4  10(eO3Cr2  4-  dood^  ) 
+  12(c12CJ2  4  d\2d\2)  +  9(021  CJ2  -f  d2\d\2  )  4  4(e3o£'i2  4  d^odvz 
+  4(co3C2i  4- do3d2i )  +  9(ci2e2i  4- d^d^i )  4- 12(c2i  021  4-  rf-jidsi 
4-  10(c3oC2i  4-  dzod2i )  4-  (C03C30  4-  do3d3o )  4  4(ci2C3o  4  d\2d^  ) 
4  10( e2ie3o  4  d2id3o)  4  20(e3oC3o  4  d3od3o)]/i 
h 

=  a — ‘  1575. 

2  •  7! 

(c)  For  b3_]  _j,  we  also  need  the  Bernstein  coefficients 


{etm}  =  (^,0.0,o) 


in  addition  to  Eq.  (18)  for  N$(t),  yielding: 

(3!)2  r _ //  7  \  2  /1\2\ 


‘4,=fK©4(;r)+4^f(rs) 


'2  7 


+  (5-5)  +  »(H-5)  +  “(5)+«(v5) 

+  - 1)  +  <5  •  I)  *  •  ■  (M)  +  “((*)’  +  (§)■) 

+  »(rS  +  8i)*(s-5)+‘(5)(5) 

+  M5-5  +  i>»(Q,  +  (8)*)]* 


3294. 


Combining  the  results  in  (a)-(c),  we  have  obtained  the  blocks  E2  and  E2  in  Eq. 
(11).  We  will  next  verify  the  correctness  of  the  blocks  D\  and  Dj . 

(d)  For  Ido  o  =  63,-  _2,  we  need  the  Bernstein  coefficients 


{Qrn}  -  (o,0,0,6),  {dpg}  -  (G<G-6-G) 


of  ( t )  in  addition  to  those  in  Eq.  (17)  of  .V;j (t),  yielding: 


'>U.«  -  '4-2  -  f  [20(g)  (g)  +  10(j)  (g)  +  4(g) 

+00+^)+i<o)+^ 


■  34S. 
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(e)  For  b3_ 2  ,  =  '.'j  _2 -  we  only  need  the  Bernstein  coefficients  { nin }  of  JV3 ( t )  in  Eq.  (20) 
as  well  as  those  of  jY3 ( t )  in  Eq.  (17),  yielding: 

. ,  . ,  3G  r 1  1 7  .  b 


b_2A-bx  .2-  7r[4  ■  6J/j  -  2.7.  -3- 

(f)  Since  the  supports  of  N3(t)  and  jV3 ( t  —  j),  for  /  >  2,  do  not  overlap,  we  have  }  ~ 
ft]  _2  =  0  for  all  >  2. 

(g)  For  fti,  o  =  6^-1*  've  need  the  Bernstein  coefficients  of  .Y3(/)  in  Eqs.  (IS)  and  (  19). 
and  those  of  A  3 ( t )  in  Eq.  (20)  and 

..  .  /4  4  2  h 


-  (o'  6’  6"  g) 


yielding: 


(3!)2  ^  A  (f  +  p)!(G-7-p)! 


'-i,o  ~  "o . - 1 


t 1  v'  v-' 

_  7!  2-J  2s  (\ 


t=0  p—0 


=  o^(2264)’ 


C!(3  -  0!/>!(3  -  />)! 


[r((~P  +  /!( d  p  +  c(Cp]h 


where  the  second  subscripts  in  and  q  in  C(m,cpq,  etc.  have  been  deleted  for  convenience, 

(h)  For  fti,  ,  =  b\  we  need  {dp?}  in  Eq.  (18),  in  Eq.  (19),  and  {c(m }.  {dpq }  in 

Eq.  (20),  yielding 


■“G-SMh-S+HMS-S) 

♦“(IXS)]* 


(i)  For  b3_ ]  2  =  b2  we  only  need  {eP9}  and  yielding: 

*„  =  lj  ,  *  .2. 

“'■2  2'“‘  7!  16  6J  2  7! 

(j)  For  j  >  3,  we  have  bi_l  =  6]  —  0,  since  the  supports  of  the  two  77-splines  .Y3 ( t ) 

and  A’3 ( f  —  j).  where  j  >  3,  do  not  overlap. 

Combining  the  results  in  (d)-(j),  we  have  obtained  the  blocks  D\  and  D[  in  Eq. 
(12)  as  well  as  the  two  zero  blocks  of  -43  /,  at  the  upper-right  and  lower-left  corners  in 

Eq.  (10).  Since  the  other  blocks  C„_  (.  D .  D  .  and  E"3  of  A  t  /,  have  been  verified  in  Ref. 
1.  we  have  now  obtained  the  formulation  of  Eq.  (10)  for  .4  ( 
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Estimation  by  Cubic  Splines 

For  noise-fme  or  very  low  noise  signals,  we  may  use  L 2  =  L2[0,d\  approximation  as 
suggested  in  Ref.  1.  Let  -  1,.  . .  ,  A\  and  0  <  T\  <  ■  ■  ■  <  <  d,  denote  the 

data  measurement,  where  the  number  N  of  sample  points  is  usually  much  larger  than 

the  number  n  of  knots  of  the  cubic  B-spline  series  Eq.  (6).  As  in  Ref.  1,  let  /(f)  denote 
the  piecewise  linear  continuous  function  on  [0,d],  linear  on  each  interval  [r,-,  r,+  i]  such 

that  /(r,)  =  /,.  This  continuous  model  of  the  discrete  signal  { r, ,  /, }  is  recommended 
only  in  a  noise-free  or  low-noise  environment.  As  in  Ref.  1  (pages  12-13),  we  consider 

the  L 2  =  Z2[0,  d]  estimation  of  /(f)  by  a  cubic  spline  curve  given  by  the  spline  series 
S3(f)  in  Eq.  (6).  Fix  the  number  n  of  knots  (n  <<  iV),  and  for  the  time  being  also  let 
h  (=  (d  —  c)/n)  be  fixed  The  L2  model  with  “smoothing  parameter”  A  considered  in 
Ref.  1  is  to  minimize  the  functional 


A'3(/i,c)=  /  /(f)-  £  CjBjit) 

JO  1  j=-2 


(It 


(22) 


■>=— 2 


over  all  c_ 2, . . . ,  cn_! ,  where 


c  = 


C-2 
C  — 1 

Cn  —  \  J 

We  will  discuss  this  model  in  details  in  the  section  entitled  Penalized  Least-Squares  Esti¬ 
mation.  Here,  the  B-splines  Bj(t)  are  given  by  Eqs.  (3)  to  (5).  Let 


n  — 1 


s3{t)  -  s3{t-,  \,h)  =  Cj Bj{t) 


i=-  2 


be  the  (unique)  minimum  solution;  that  is, 

E3{h)  :=  I\3(h,c)  =  min  J\3(/i, c) 


(23) 


(24) 


where  c  =  [c_2,  •  •  •  ,c„_i]r  depends  on  A  and  h.  Then  c  =  c(A,  h)  can  be  determined  In- 
solving  the  normal  equations 


n  —  1 


f>ztjc}  +  Xc,  =  /h°  1  =  -2 . u-\ 

j=- 2 


( 23 1 


where 


n,  =  f 

./n 


f(t)Bx(t)df. 


20  \ 
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In  matrix  form.  Eq.  (24)  becomes 


(Ai,/i  +  A/b+2)c  —  (27) 

where  A3j,  is  the  coefficient  matrix  given  by  Eq.  (10).  I„+2  is  the  (n+2)  x(» +2)  identity 
matrix,  and 


tt  = 


fl- 2 


If, 


h,n—  1  J 


is  the  (n  +  2)-dimensional  data  vector.  To  compute  f£,  a  change  of  variables  in  Eq.  (25) 
yields: 


and  for  i  =  0, . . .  ,  n  —  1, 


fl- 2  =  h 

[  f{h(t  -  n)  +  d)X(t)dt, 

Jo 

(28) 

fl-y  =  h 

[  f(t>(t  -  n)  +  d)X(t)dt, 

Jo 

(29) 

fl,  =  h  f 

Jo 

-t 

f(h(t  -  n  +  i)  4-  d)X(t)dt. 

(30) 

Note  that  in  Eq.  (30)  the  integral  is  taken  over  the  interval  [0,4]  for  i  =  0 . n  —  4, 

but  is  taken  over  a  smaller  interval  for  i  =  n  —  3,  n  —  2,  n  -  1.  According  to  page  32  and 
Theorem  4  on  page  25  in  Ref.  1,  to  determine  the  signal  onset  to,  we  must  determine 
the  set  H  =  {/?}  such  that  each  h  satisfies: 


E3(h)  =  min  E3(h)  (31) 

h  ^0 

where  E3(h)  —  I\3(h,c(X,  h))  is  defined  in  Eq.  (24).  Here,  since  only  noise-free  or  low- 
noise  environment  is  considered,  A  is  fixed  and  in  fact  may  be  set  to  be  zero  for  noise- 
free  signals  or  a  very  small  positive  number  for  signals  with  low-noise  contamination. 
Then  the  signal  onset  (or  time-of-arrival)  is  given  by 


to  =  d  —  nh*  (32) 

where  h*  is  the  smallest  number  in  H  =  {^}. 

In  Ref.  1  (pages  25-20),  to  determine  the  set.  H  =  {h}.  we  set  D  =  (1  /h)A3  i,  and 
b  =  (l//))f°;  i.e.. 


\D  = 

[— 63  ] 

h  ,}\ 

lb  = 

— 2<i.ji<n  —  1 

~2<t<r)-l 

[which  is  essentially  Eq.  (74)  in  Ref.  1],  and  write 

B  =  PAP1 


(331 


341 
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where 


A  = 


■M 


o 


(35) 


O  -^n  +  2 

with  Aj  >  •  •  ■  >  An+2  >  0  being  the  eigenvalues  of  B  and 

P  =  [uj  . . .  un+2]  (36) 

with  u,  being  the  eigenvector  of  B  corresponding  to  A,  and  having  unit  length.  Also,  set 

o 


r  = 


hXi*f  X 


o 


h\n+2 + A 


(37) 


and 


b  =  Prb  = 


bx 


bn 


+2  J 


(38) 


where  b  is  defined  in  Eq.  (33).  Then  Ez(h)  in  Eq.  (24)  becomes 

rd  «+2  ,2 


fa  1Z  h 2 

E^=L 


(39) 


where  b,  is  the  i-th  component  of  b  in  Eq.  (38).  Hence,  each  h  is  obtained  by  finding  an 
absolute  maximum  of  the  expression 


n+2 

™  =  E  nd 


i=  1 


h\i  "f*  A 


(40) 


where 


d,  =  di(h)  =  /i6,.  (41) 

We  summarize  the  above  discussion  as  follows. 

Algorithm  I  (Cubic  Spline  Estimation  of  Signal  Onset  in  Noise- Free  or  Low- 

Noise  Situations) 

(1°)  Choose  A  =  0  (if  the  signal  is  noise-free)  or  a  positive  but  small  value  of  A  (for  low- 
noise  signal).  (In  the  next  two  sections,  we  will  discuss  procedures  for  estimating  A 
when  the  signal  is  fairly  noisy.)  Also,  choose  a  positive  integer  n  such  as  tj  =  10  or 
anything  larger.  (The  dimension  of  the  matrix  will  be  (n  +  2)  by  (n  +  2).) 

(2°)  Compute  the  eigenvalue-eigenvector  pairs  (A,,u,),  ?  =  1 _ _ n  +  2.  of  the  matrix 

B  in  Eq.  (33)  using  the  entries  bf  of  A^h  given  by  Eqs.  (10)  to  (15),  where  u,  is 
normalized  to  have  unit  length.  Note  that  the  variable  h  does  not  appear  at  this 
stage. 
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(3°)  Form  the  matrices  P  and  F. 

(4°)  Compute  the  data  vector  f°  =  [/12,-  •  •  '/n-i]T  hi  Eqs.  (28)  to  (30)  and 


■  dx  ■ 

■/“  2- 

d  = 

=  PT 

.  hn+ 2  . 

1- 

(5°)  Plot  the  curve  for  T3(h)  in  Eq.  (40)  and  determine  its  set  of  absolute  maxima 

H  =  {h}. 

(6°)  Determine  the  smallest  value  h*  in  H.  Then  the  signal  onset  is  t0  =  d  —  nh* . 
(7°)  Compute 


c  =  /Tb  (44) 

where  c  =  [c_2,  •  • . .  en-i  ]  r • 

(8°)  Plot  the  cubic  spline  curve 

71—1 

s3(t)  =  ^2  c,B,(t)  (45) 

i=-  2 

where 

=  N3  (f  ~  <o)) 

and 

Bj(t)  =  iV3  (t  -  t0)  - 

for  j  =  0, . . .  ,n  —  2.  [See  Eqs.  (3)  to  (5)].  Then  S3(t)  is  the  cubic  spline  estimate  of 
the  underwater  acoustic  signal  with  signal  onset  at  t0. 

PENALIZED  LEAST-SQUARES  ESTIMATION 

We  are  interested  in  fitting  a  set  of  noisy  data  {r,, /,}.  i  =  1 . AT,  where 

0  <  fi  <  •  •  •  <  7-y  <  d,  by  a  “smooth”  curve  g(t )  on  [c.  d],  which  is  continuous  on 
(0,d|  such  that  g(t)  =  0  for  0  <  t  <  c  but  is  not  identically  zero  on  [0,c'[  for  any  c'  >  c. 
Hence,  c  =  to  represents  the  signal  onset  of  the  underwater  acoustic  signal.  We  assume 
that  the  actual  signal  is  given  by  a  function  /*(/)  in  the  sense  that 

fi  =  f*(Ti)  +  c,  (4G) 

where  c,’s  are  uncorrelated  noise  processes  with  zero  mean  and  positive  variances:  that 
is. 


S(£i)  =  0  and  £{s[  £,)  = 
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where  a,  >  0 ,6tJ  is  the  Kronecker  delta  (which  is  defined  to  be  1  for  i  =  j  and  0  for 
i  ^  j),  and  £  is  the  expectation  operator.  Our  choice  of  the  “smoothing"  data-fitting 
curve  g(t)  is  motivated  by  considering  nonparametric  regressions  in  the  sense  that  no 
mathematical  model  of  the  smooth  function  g(t)  is  a  priori  assumed.  Here,  "smooth¬ 
ness"  means  that  g(t)  belongs  to  the  class 

W™ o  =  {g  £  Cm_1[c,d]:  g(c)  =  0,</m)(t)  square  integrable}  (47) 

of  functions  on  [c,  </],  having  zero  extension  to  the  whole  interval  [0,  d],  where  0  <  c  =  to- 
The  subscripts  2  and  0  indicate  that  least-squares  measurement  will  be  used  and  the 
functions  all  vanish  at  c  —  tQ ,  respectively,  and  m  is  any  preassigned  non-negative  in¬ 
teger  of  our  choice  of  the  order  smoothness.  Of  course,  Cm~l{c,<J\  denotes,  as  usual,  the 

collection  of  all  functions  g(t)  such  that  g(t),  g'{t)-,  ■  ■  ■ ,  *(t)  arc  all  continuous  on 

[c.d\.  Note  that  ir2"'0  is  a  proper  subspace  of  the  so-called  Sobolev  space  W™  on  [c,  d] 
where  the  condition  g(c )  =  0  is  not  assumed.  It  must  be  emphasized  that  any  function 
g (t)  in  W™o  is  extended  to  [0,d]  by  setting  g(t)  =  0  for  0  <  t  <  c  to  represent  the  under¬ 
water  acoustic  signal.  (We  remark  that  our  algorithms  in  this  report  can  be  modified  to 
study  parametric  regressions  using  sinusoidal  “wave”  functions  although  the  calculations 
would  have  to  become  much  more  complicated.) 

Since  the  data  (r,-,/,)  are  noisy,  no  reliable  way  is  available  to  model  the  data  by  a 

continuous  function  /(f)  as  before,  so  that  the  L 2  norm  cannot  be  used.  Instead,  we  will 
use  the  C2(w)  norm  defined  by 


llslb  =  llg(U2(w)  = 


(48) 


where  each  iC{  is  strictly  positive  with  uq  +  •  •  •  +  wn  =  N  and  g  =  {?«}•  ^  9 >  =  9(Ti) 
for  some  function  g(t),  we  will  simply  use  the  notation  ||^||2  =  ||g||2-  The  sequence  w  = 
{»’,}  is  called  the  weight  of  the  i 2  =  ^2(wl  norm.  When  the  variances  cr2’s  are  all  non¬ 
zero  as  we  have  assumed  here,  the  optimal  weight  is 


Ncr 


-2 


wi  =  — T5 
al 


-f - her 


-2  ’ 
N 


(49) 


the  reciprocals  of  the  variances  so  normalized  that  uq  +  ■  ■  ■  +  u\\  =  N.  (See  [7],  page  18). 
Hence,  if  the  variances  can  be  measured  and  each  a 2  is  nonzero  we  will  always  choose  w, 
as  in  Eq.  (49);  but  if  they  cannot  be  estimated,  we  will  simply  set  w,  =  1  for  all  i.  In 
the  case  when  7V0 , 1  <  JV0  <  N,  of  the  <72’s  are  zero,  then  the  weights  corresponding  to 
the  zero  erf' s  are  all  equal  to  iV./V0'’1,  while  the  weight  0  is  assigned  to  the  remaining  ones 
corresponding  to  nonzero  af's.  This  weight  assignment  is  simply  a  consequence  of  the 
limiting  process  from  Eq.  (49). 

Given  a  data  sequence  {r,,/,},  we  are  interested  in  studying  the  minimization  of  the 
functional: 


C-r,(<7,  A)  =  ||9  -  f|||  +  A  f 

i  y  r1 

=  Jr  ^2(g{T,)  -  /,)*«',  +  a  /  ( g(m'tt)V<lt 

_ 1  J  C 


(50) 
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where  A  >  0  is  called  the  smoothing  parameter.  Note  that  A  —  0  corresponds  to  strictly 
least-squares  data  fitting,  but  for  very  large  values  of  A,  the  term  jj g  -  /l|2  becomes  neg¬ 
ligible.  so  that  the  functional  Ghlg,  A)  is  dominated  by  the  “smoothing"  operator  defined 
by  the  square  of  the  ?7?- 1 h  order  derivative.  Our  objective  is  to  find  the  best  data  fitting 
curve  while  guaranteeing  certain  degree  of  smoothness  governed  by  in  and  A:  that  is.  we 
are  interested  in  obtaining  a  curve  g(t)  that  minimizes  the  functional  Gft(g.  A),  or  some 
modified  formulations  of  it,  over  all  functions  g(t)  in  U  "*0  and  all  smoothing  parameters 

A. 

Since  functions  from  H’2"o  are  used,  they  have  continuous  extensions  from  [e.  c/j  to 
[M  such  that  the  extended  values  are  identically  zero.  Hence,  an  “extremal  function" 
in  the  above  discussion  represents  the  acoustic  signal  with  signal  onset  at  c  —  t0 .  Conse¬ 
quently.  determining  the  set  H  =  {/i}  such  that  each  h  >  0  satisfies 

min  G-k  (o,  A* )  =  min  min  Gh(q,\*),  (51) 

g  €  IV™0  ’  h>0 

where  A*  is  characterized  by  minimization  of  the  weighted  mean-squares  error  in  approx¬ 
imating  the  given  data  (or  some  modification  of  it)  yields  the  signal  onset  t0  =  d  —  nli *. 

where  h *  is  the  minimum  value  among  all  h  in  H  (See  [1],  pages  25-2C.  See  subsection 
entitled  The  Generalized  Cross-Validation  Function  in  this  report  for  a  more  precise'  no¬ 
tion  of  A* ). 

Spline  Solution  of  the  Minimization  Problem 

In  the  following  we  will  motivate  our  approach  to  the  study  of  our  extremal  prob¬ 
lems  using  splines  by  studying  the  minimization  problem: 

Em(h,X)  =  min  Gh(g,X)  (52) 

9^  VV2,0 

for  any  fixed  positive  values  of  h  and  A.  Here  the  data  {r,  ,  f, }  are  to  be  estimated.  We 
first  remark  that  it  is  a  standard  mathematical  argument  (using  the  completeness  of 
H-)  that  a  function  g*(t)  in  IV2m0  exists,  such  that 

Gh{g*,  A)  =  Em(h,  A).  (53) 

Next,  we  will  see  that  any  extremal  function  g*{t)  [satisfying  Eq.  (53)]  must  also  satisfy 
the  condition: 


i: 


g*{m\t)C{m)(t)dt  =  0,  all  £{t)  €  V2m0, 


(54! 


where  V2mn  is  defined  by 

vVTo  =  {^  G  w2y.  C(rt)  =  0.  /  =  1 . V}. 

To  verify  this  fact,  let  us  consider  the  new  data 

=  <7 *(74  ) 

and  the  corresponding  subcollection 

=  iff  €  H'”’:  g(rt)  =  ct.  f  =  1 . V  } 


(  .jo 
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of  functions  in  W^o  that  interpolate  the  same  data  as  g*{t).  Hence,  we  have 

0<Gk(g,\)-Gk(g*,\) 


for  all  g  €  particular,  for  any  £  £  \ 'jm0  and  any  real  number  s.  since  g*  +  si  is  in 


W2m0,  the  function 


rd 

F(s)  =  j  (g^m\t)  +  ^m\t))2dt, 


which  is  a  real-valued  differentiable  function  in  s,  has  a  relative  minimum  at  .$  —  0.  so 
that  F'(0)  =  0.  In  other  words,  we  have  verified  Eq.  (54). 

Next,  we  will  verify  that  g*(t)  is  essentially  unique  if  the  number  N  of  r,’s  that  lie 
in  the  interval  [c,  d]  is  at  least  m,  in  the  sense  that  for  any  g(t)  in  IV that  also  satisfies 
Eq.  (54)  as  g*{t)  does,  we  have  g(t)  =  g*(t).  Indeed,  if  g(t)  is  in  W™c ,  then  g{t)  —  y*(t)  is 

in  V2m0  and  can  be  used  as  a  function  £(£)  in  Eq.  (54).  Hence,  since  both  g(t )  and  g*(t) 
satisfy  Eq.  (54),  their  difference  also  satisfies  Eq.  (54),  so  that 


rd 

0  =  J  (s<m,(t)  - 


which  implies  that  g(m\t)  —  g*(m\t)  =  0,  or  g(t)  —  g*(t)  is  a  polynomial  with  degree  at 

most  m  on  [c,  d].  But  since  this  polynomial  has  N  zeros  and  N  >  m,  we  may  conclude 
that  g(t)  —  g*{t)  =  0. 

The  above  observation  allows  us  to  characterize  g*(t)  as  follows.  Let  S(t)  be  a  func¬ 
tion  in  C2m~2\c,  dl  whose  restriction  to  each  of  the  intervals  [c,  r..  ~  l,[r.  ~ 

r^v— ~+2], ....  [r/v_i,r/v],  and  [r^/,d]  is  a  polynomial  with  degree  at  most  2m  — 1;  i.e..  S(t.) 
is  a  spline  of  degree  2m  —  1  in  [c,  dj.  Let  us  first  assume  that  c  ^  rAf_^  +  1  and  d  ^  r.v,  so 
that  there  aie  precisely  N  interior  knots:  rN_~  . . . ,  in  [c,d],  and  this  implies  that 
S{t)  has  2m  +  N  free  parameters,  2m  from  the  polynomial  piece  on  [c,  r/v_~+1],  and  1 

from  each  of  the  N  interior  knots  in  [c,  e/j .  (The  argument  will  be  similar  if  one  or  both 
of  and  r,v  should  become  non-interior.)  Hence,  we  may  construct  the  (unique) 

spline  S(t)  that  satisfies  the  interpolatory  conditions: 


(i)  5(c)  =  0, 

(ii)  5(7-,)  =  z,.  for  i  =  N  -  N  +  1 . JV. 

(hi)  5(m)(c)  =  •  •  •  =  5(2m_2)(c)  =  0,  and 
(iv)  S{m)(d)  =  ■  •  ■  =  5(2m -0((/)  =  0; 


a  total  of  1  4-  A’  4-  (m  —  1)  +  m  =  2m  4-  .V  conditoins.  If  rv_-  =  c.  then  the  condition 
(i)  is  not  necessary  since  it  is  included  in  (ii);  and  if  r,\  —  d.  then  it  will  be  seen  latei 
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that  the  condition  S(2r"  ])(d)  =  0  in  (iv)  can  be  deleted.  In  any  case,  the  spline  S(t) 
satisfying  (i)-(iv)  (or  the  corresponding  ones  if  r,\-  =  d)  exists  and  is  actually  unique. 

By  setting  S[t)  —  0  for  0  <  /  <  c,  it  follows  from  condition  (i)  that  S(t)  is  in  ITT',,. 

(Actually  it  is  somehwat  smoother,  being  in  C2m~2[c.  d],  since  2m  —  2  >  in  —  1  for  ///  >  1.) 
We  will  now  verify  that,  in  fact,  g*(t)  =  S(i)  by  showing  that  S(t)  satisfies  Eq.  ( 34 )  as 

g*(t)  does.  [Observe  that  S(t)  is  in  I V’2"'0  by  the  mterpolatory  condition  (ii).j  We  remark 

that  S(t)  is  not  a  natural  spline,  since  the  condition  *(r)  =  0  must  lx-  dropped 

and  an  interpolation  condition  S(d)  =  +  i  is  required  for  it  to  qualify  as  a  natural 

spline  of  order  2 m. 

To  verify  Eq.  (54).  let  (( t )  be  any  function  in  E/",.  Then  in  addition  to  the  interpo¬ 
lator}'  conditions  f(r, )  =  0,  we  also  have  f(c)  =  0.  For  notational  convenience,  we  set 
c  =  .r„,rv_~  +  1  =  j't - -  r,v  =  x~,  and  d  =  x~+r  so  that 

f(.r,)  =  0,  i  =  0 . X  (5S) 

and 

S(.r,)  =  r,,  i  =  0 . -V. 

,  S(ml(.r„)  =  •••  =  =  0.  (59) 

U("l<W  =  -''  =  S'2'"-|’<Jv+I>  =  a 

By  applying  integration  by  parts,  we  then  have,  for  any  ((t)  in  V.™Q: 


cn  -v+i 

/  S(m)(t)f.{m)(t)dt  =  Y 

Jc  )  =  1  d-r.-i 


Sim)(t)(im)(t)dt 


=  { S( \.r  ,)  —  S( ”' * ( .r i - 1  )f ( ~ 1 1  ( x, _  j ) 

i 

-  f  '  S(m+1){t)C{m~i)(t)dt}  =  ■■■ 


•V  -f  1  ,  m  —  1 

-  Y  |  ^(-lH[5(”,+-')(.r,K<m-^,)(.rI)  -  5(m+j)(x.--iK(n,">-lVi-i 

i=i  t  J=0 

+  (-l  )m  '  Si2",](t)C(t)dt  j 


+  (-!)' 


^(-lF[S(n,+2)(.r~  +  iK,"|-2-,)(.r~+i)-S,"‘+2,(.r()K,"'-2-1)(.r())] 


'V+'  rr, 


+  (-!)' 


"  T  1  [r< 

Y  /  SV2m\t)f(t)dt  =  0. 

,=i 


by  telescoping,  applying  E([s.  (5S)  and  (59).  and  using  the  fact  that  S{2n,)(t)  —  0  on 
each  _ | .  ./*,).  respectively.  This  complets  the  proof  of  the  claim  that  </*(t)  ~  S{t ). 
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Modification  of  the  Extremal  Problem 

We  have  just  proved  that  spline  functions  with  knots  given  by  the  sample 
points  { rt }  where  the  data  {/,}  are  taken  must  be  used  to  yield  optimal  penalized 
estimation  of  the  noisy  data  {r,,/,}  in  nonparametric  regression.  However,  since  the 
points  {  t,  }  are  chosen  without  any  knowledge  of  the  signal  onset  t o,  they  cannot  be 
used  as  knots  to  estimate  t0.  Note  that  t0  is  a  variable  in  our  mathematical  model. 

In  addition,  the  number  N  of  measurements  is  usually  much  larger  than  the  realistic 
number  n  of  knots.  In  order  not  to  get  too  far  away  from  nonparametric  optimal  pe¬ 
nalized  estimation,  we  therefore  restrict  our  attention  to  the  subspace  0  with  basis 

{Bk'th.-k+ . . Bk  th  „-i }  of  o  where  k  =  2m  —  1.  That  is.  every  spline  subspace 

function  in  0  is  a  spline  series 


Sk(t) 


n  —  1 


yi  cjBk, 

>=-*+ 1 


(60) 


at  discussed  before.  Here, 


t/j!  —  •  —  ^—1  —  (^1) 

•  k 

with  0  <  t0  =  c,  <  d  and  tn  =  d ,  is  the  knot  sequence  of  the  space  St;>  0,  as  discussed 
in  the  Introduction.  Observe  the  effect  that  if  n  is  chosen  to  be  N  +  1  and  the  r,'s  arc 
also  equally  spaced,  then  the  best  choice  of  t0  will  enable  us  to  reproduce  the  nonpara¬ 
metric  optimal  penalized  estimation.  It  must  be  emphasized  that  only  odd  degree  splines 
should  be  vised  to  retain  the  spirit  of  nonparametric  optimization;  hence,  corresponding 
to  m  =  1,2,  degrees  k  =  1  and  3  (i.e.  linear  and  cubic  splines)  will  be  most  useful.  For 
g(t )  =  Sk(t)  given  by  the  spline  series  in  Eq.  (60),  the  functional  Gh(g ,  A)  in  Eq.  (50) 
that  we  are  interested  in  minimizing  becomes 


Gk{Sk,  A)  =  US,  -  f||i  +  A  J\s<km](t))2dt 


C 

1  N  n-1  2 

=  jy  ~  U'1 


i=l  L  ;=-*+ 1 

n  — 1  n  —  1 


+  *  E  E  c>c>  f 

Lit-  L  «  l  Jc 


*=  — fc+1  j  =  -k+ 


(62) 


The  only  difference  of  this  from  the  minimization  problem  Eq.  (52)  is  that  we  now  mini¬ 
mize  over  all  Sk(t)  in  the  subspace  0  instead  of  over  all  g{t)  the  entire  space  in 

other  words,  the  minimization  is  over  all  coefficients  c_,+  i . c„_).  Hence,  as  in  Ref.  1 

(pages  9-10),  to  find  the  spline  series 


n  —  1 

sk(f)=  c]B^tk.Af) 

that  satisfies 

Gh(S;,  A)  =  Ek(h.  A). 


(63) 


(64) 
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where 


Ek{h,X)=  min  Gh(Sk.  A), 

S‘eSt\.o 

(Go) 

we  simply  differentiate  Gh(Sk,X)  in  Eq.  (G2)  with  respect  to  each  of  o_t+1..  . 
yield  the  normal  equation.*: 

.  .  ,C„_1  to 

where 

11-1 

y.  (°k.,j  +  NXb™tJ)c*  =  /jt.i,  i  =  -k-  +  1 . v  —  l. 

j  —  —  1 

(GC ) 

N 

ak.,j  =  y,Bk,th,i(Tt)Bk.tH.j(T()u't , 
e=\ 

( G7) 

h,,  =  [ 

(GS) 

and 

N 

h,i  =  y  ,i(  T()w(- 

(= i 

(09) 

By  setting 

Ak.h  =  —  t+l  <i.J<n  —  1 

(70) 

Bk.h  =  [&fc,i>]-fc+1<i.j<r«-l 

(71) 

c-*+i 

.fh-k+\ 

c*  — 

f h  = 

.  c»-i  j 

J h ,  n  —  1 

the  matrix  formulation  of  the  normal  equations  in  Eq.  (GO)  becomes: 

(Ak.h  +  jVAZ?°Jc*  -  f,(.  (73) 

In  computation,  we  first  determine  the  N  x  (n  +  k  —  1)  matrix 

Bk.h  =  [Bk,th  ,>(r,)]  i  < ,  <  \  •  (74) 

-/t  +  l<7<n-  I 

Then  setting 


"4 

o1 

■  f,  ■ 

w  = 

f  = 

.0 

fC  V  _ 

.  /  V  . 

IS 
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where  {/,}  is  the  noisy  data  taken  at  {r,},  the  normal  equations  in  Eq.  (G6),  or  equiva¬ 
lently  Eq.  (73),  become 

(Bl,WBk,k  +  (70) 

It  should  be  remarked  that  the  (n  +  k  -  1)  x  (n  -f-  k  ~  1 )  matrix  B°k  h  can  be  pre-computed 
as  what  we  will  do  in  the  next  sub-section.  There,  we  will  see  that  in  general  this  matrix 
Bk  h  is  fairly  complicated.  Since  it  has  to  be  multiplied  by  the  smoothing  parameter  A 
which  still  has  to  be  estimated,  inverting  the  coefficient  matrix  of  the  unknown  vector 
c*  in  Eq.  (76)  to  find  c*  causes  some  complication.  For  this  reason,  we  proposed  [Eq. 

(27)  on  pagelO  of  Ref.  1]  to  replace  the  matrix  Bkh  by  the  identity  matrix  I;  that  is.  the 
modification 

(Blj,WBk,h  +  N  AJ)c*  =  Bj  b\Yt  (77) 

of  Eq.  (7G)  will  also  be  considered. 

Computation  of  the  Matrix  Bk  k 

If  we  wish  to  use  the  original  normal  equations  in  Eq.  (66),  or  equivalently  (76),  we 
must  pre-compute  the  matrix  B°k  h.  For  linear  splines  (i.e  m  -  k  =  1),  this  can  be  done 
easily.  Indeed,  by  Eq.  (2)  for  k  =  1,  we  have  B\M  }{t)  =  ( l/h)N[[(l/h)(t  -  c)  -  j]  where 


N'(t) 


1  for  0  <  t  <  1 

-1  for  1  <  t  <  2, 

0  otherwise 


so  that 


b 


i  M  ~ 


l 

1 

h 


d 


B\ 

f 

N'(t  -  i)N'(t 

-  j  )dt 

Jo 

2 

h 

for 

i  =  >1 

1 

h 

for 

*  =  j , 

1 

~~R 

for 

1*  -  J 1 

0 

otherwise. 

0  <  i  <  n  —  2, 
i  —  n  —  1, 


That  is.  the  matrix  B°  h  in  Eq.  (76)  for  k  =  1  is  given  by 


B°  -  - 

u\,h  -  . 


2 

-1 

o 


-1 

o 


-1 


o 

2  -1 

-1  1 


(78) 


(79 
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For  cubic  splines  (i.e.  k  =  3,  m  =  2),  the  computation  is  much  more  involved.  We 
will  again  use  Bernstein  representations.  Hence,  from  Figs.  2  and  3,  and  the  correspond¬ 
ing  expressions  that  follow  the  figures,  the  Bernstein  representations  of  N"(t ),  X^t ). 
and  N"(t)  are  piecewise  linear  polynomials  given  by  the  Bernstein  representations  in 
Fig.  4  below: 


Fig.  4  -  Bernstein  representations  of  and 

Hence,  by  applying  the  integration  formula  in  Ecp  (1C),  we  have  the  following  re¬ 
sults: 

(a)  For  i  =  j  =  —2,  we  have 

_  24  288 

“  h*  “  2  x  3 lh3' 

(b)  For  i  =  —2,;  =  —  1  or  i  =  —  l,j  =  —2, 

1  f2  ~ 

*3;— 2,-1  =  ^3 ;  —  1 ,  —  2  =  ^3  J  N"(t.)N”(t)dt 

27  _  81 

~  ~4 h*  ~  2  x  3!/W 

(c)  For  /  =  —  2. j  =  0  or  ?  =  0<v/  =  —2.  we  have 

1  f2  ~ 

^3 ;  —  2 . 0  =  ^3;0.-2  =  pj  /  -\”(/ )  A3  (  t  )(It 

1  _  12 

“  _/C  -  ~  2  X  3W!  ' 

(d)  For  7  =  -2.j  =  1  or  7  =  l,j  =  -2.  we  have 
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(r)  For  i 

(f)  For  i 

(g)  For  i 

(h)  For  i 

(i)  For  i 

(j )  For  ?  = 

(k)  For  ?  = 

(l)  For  >  = 


1  f  ~ 

h  -2,1  =  *3;1  .-2  =  N3(t)N3{t  ~  1  )dt 

1  3 

4 h3  ~  2  x  3\h3  ‘ 

_2,j  >  2,  or  i  >2 ,j  =  —2,  we  have 

^3;-2j  =  &3;i,-2  =  0. 

— 1,7  =  —1,  we  have 


9  54 


2h3  2  x  3!/i3  ' 

“1,7  —  9  or  i  =  0,  j  =  —  1,  we  have 

1  r3  ~ 

63;-l,0  =  63;0,-l  =  ^3  J 

4  16 

3/r3  ~  ~2x  3!/i3  ' 

-1,7  =  1  or  i  =  1,  j  =  — 1,  we  have 

&3;-l,l  =^3;l,-l  = 

1  1 
12/i3  ~"2x  3!/i3  ’ 

-1,7  =  2  or  i  =  2,7  =  -1,  we  have 

63;-l,2  =  &3;2,-l  =  ^3  [  ^3  (  W V  ~  2)dt 

1  _  2 
~  6/j3  ~  2  x  3!/?3 ' 

—  1,7  >  3  or  i  >  3,j  =  —  1,  we  have 


^3;- 1.7  =  &3;., -I  =  0. 
7  and  /  =  0, . . . ,  n  —  4  we  have 


32 


2  x  3!/i3 


7  =  n  -  3,  we  have 
b\ 


3 ;  ri  —  3 , »» —  3  — 


Fjf 


2S 


3/i 3  2x3!/f» 
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(m)  For  i  = 


(n)  For  i  = 


(o)  For  i  = 


(p)  For  i  = 


(q)  For  i  = 


(r)  For  /  — 


( s )  For  / 


( t )  For  /  - 


<  u )  For  |/  - 
Hem<\ 
unit  rix 


j  =  n  —  2,  we  have 


&3;ii-2.n-2  - 
/  =7?  —  l.  we  have 


1G 


3/; 3  2  x  3!/i  * ' 


4 


2  x  3!// 3 ' 

0 _ _  7?  —  4,  j  =  7  +  1,  or  ?'  =  j  +  1.  j  —  0, . . . ,  77 —  4.  we  have 


»+]  =  °3-,j+  t 


-  1  )M 


is 


2/i 3  2x3  !/i3' 

77-3  .J  =  7  +  1  =  77  -  2.  Or  7  =77  —  2.  j  =  77  —  1,  WO  llllVe 

1 

^3;  n  —3 , n  —  2  =  &3:n-2,n-3  ~  J  3  (  ^  ) -^  3  (  ^  —  1  )f^ 


1 


12 


P  2  x  3 !  /7 3 

7)  -  2.  j  =  7  +  1  =  77  —  1,  or  7  =  ?7  -  , .  /  —  77  -  2,  we  have 

c2 


1  [ 

3;  ti  —2 , .1—  i  —  ^3;  n  —  1 ,  n  —  2  ~  J  ^ 


"'(t)NZ(t-  \)dt 


2/73  2x3!/t3' 

,77  —  4  and  j  =7+2:  or  7  =  j  +  2  and  j  =  0, . . . ,  ??  —  4,  we  have 

(2  -  2  +  2  -  2)  =  0. 


G/i3 


—  3  and  j  =  77  —  1,  or  7  =  72  —  1  and  j  =  n  —  3,  we  have 


K.j  =  p  J  NSmUt  -  2W»  =  ^7<2 -21-0. 

0 . 77  —  4  and  j  =  ?  +  3.  or  7  =  j  4  3  and  j  =  0 . 77  —  4.  we  have 

“  F  if  =  cP  = 


-  /|  >  4.  we  h;vve  6t;l7  =  0. 

summarizing  the  results  in  (a)  (u).  we  have  obtained  the  (77  +  2)  l>v 
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where 


3,h 


2  x  3!/i3 


\b°n 

i>?2 

0  ' 

b°2X 

^22 

^23 

(80) 

0 

^32 

b°33j 

b 


o 

11 


b 


o 

22 


288 

-81 

-81 

54 

> 

1.0 

b  12  — 

-12 

-16  - 

'  32 

-18 

0 

2 

-1 

-18 

• 

• 

0 

0  ' 

• 

• 

•  • .  *2 

2  * 

• 

• 

. 

'  • .  '  0 

0 

• 

• 

•  -18 

3  0 


2  0  -18  32 J 


0  ...  O' 

0  ...  0 


b°23 


2 

0 


L  —18 


O  1 

0  0 

2  0  ’ 
0  2. 


28 

-12 

0 


-12  0 
16  -6 
-6  4 


are  2  x  2, 2  x  (n  -  3),  (n  —  3)  x  (n —  3),  (n  —  3)  x  3,  and  3x3  blocks  and  621  =  ^12  1  ^32  =  ^23  ■ 
In  general,  to  compute  Bk  h  for  k  =  5,  7, ,  one  has  first  to  determine  the  corre- 


*  rv  ,  . 

sponding  £?-splines  Nkj(t)  with  £  stacked  knots  at  the  origin,  £  =  2, . . . ,  k  —  1  [where 
Nki(t)  =  Nk(t)],  then  compute  their  mth  derivatives  where  m  =  ( k  +  l)/2,  and  finally 
evaluate  the  integrals: 

■1  />k- f  1 

h,,  =  jf  (si) 

for  —k  +  1  <  ?,  j  <  —1; 

=  (82) 

for  —  fc  +  l<i<— 1  and  0  <  j  <  n  —  1.  or  0  <  i  <  r? .  —  1  and  —  A-  +  1  <  j  <  —1:  and 

bk''J  =  ^l  N{™){t-l)N*m){t-j)<It  (S3) 


for  0  <  <  n  —  1.  Of  course,  the  formula  in  Eq.  (16)  should  be  used  to  facilitate  the 

integration  procedures. 


The  Generalized  Cross-Validation  Function 


Let  us  now  return  to  the  normal  equations  in  Eq.  (76)  and  the  modification  in 
Eq.  (77)  with  Bkh  replaced  by  he  identity  matrix  I.  Note  that  the  matrix  B^j, 
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[■Ofc,t*.j(7‘t)]  is  simple  to  evaluate  and  the  weight,  matrix  IV  and  data  vector  f  are  already 
predetermined  and  given.  Hence,  it  is  simple  to  determine  c*  =  c*(/>.  A)  for  any  fixed 
values  of  h  and  A  by  using  the  Moore-Penrose  pseudoinverse  in  Ref.  1  (pages  1G-20).  We 
will  always  fix  h  >  0  in  the  following  discussion.  The  objective  is  to  determine  an  opti¬ 
mal  value  of  the  smoothing  parameter  A  in  terms  of  the  noise  variances  of  the  given  data 
f  =  {/,}.  Let  c*  =  c*(A)  =  [clit+1  . .  .c*_,]7  be  the  solution  of  Eq.  (7G)orEq.  (77) 
using  the  Moore-Penrose  pseudoinverse 


M(X)  =  l(BihW^h+NXD^)+  forE(l(7G)  (8- 

l  {Blh\VBk  h+.WI)+  for  Eq.  (77) 

in  the  sense  that  c* ( A )  =  A/(A)J9^  ^IVf  is  the  unique  solution  of  Eq.  (7G)  or  (77)  with 
minimum  £2-norm.  (See  page  16  in  Ref.  1).  In  other  words. 


S*{U  A)=  5]  c*j(\)Bk<thJ(t) 


J=-k+i 


is  the  (optimal)  penalized.  least-squares  estimator  of  the  noisy  data  f  =  {/,}  taken  at  the 
sample  nodes  {r,}.  To  estimate  the  optimal  smoothing  parameter  A,  we  investigate  how 
well  S*(t\  A 't  pproximates  the  actual  (unknown)  signal  {/*(r,)}  [when  the  noise  {c,}  is 
removed  from  the  measured  data  {/,};  see  Eq.  (46)  for  the  definition  of  /*(r,)]:  that  is, 
we  are  interested  in  the  size  of  the  quantity 


1  ‘ 

rf(A)  =  -J(q(5V,;A)-r(r,))2. 


Since  Tf(\)  depends  on  the  noisy  data  f,  we  must,  in  fact,  investigate  its  mathematical 
expectation  £7f(A).  E  seems  then  very  reasonable  to  characterize  the  optimal  smoothing 
parameter  A*  =  A^-  as  one  that  minimizes  ETf(A),  namely: 

<fTf(A*v)  =  min<f7f(A).  (87) 

A>0 

However,  this  is  an  impossible  task  since  the  actual  signal  {/*(*,)}  in  the  definition  of 
Tf ( A )  is  unknown.  For  this  reason,  the  alternative  quantity 


1  N 

Tf{\)  =  -y  U’*(5’*(r<:  -M  -  />  )' 


which  measure  the  error  in  approximating  the  dat  a  r  -  {/,}  by  fbe  optimal  penalized 
least-square*  estimator  5*(t;  A)  must  be  taken  into  consideration.  In  fact,  it  will  be  seen 

that  a  weighted  quantity  of  Tf(A)  given  by 


Vv(A) = 


7f(A) 

(1  -  I.v(A))2* 


where 


L  v ( A )  =  -  Tvm-c[Dk,hM(\)D[h]V] 
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with  M( A)  defined  in  Eq.  (84),  shall  give  some  information  on  the  optimal  A*v.  The  util¬ 
ity  of  the  weight  (1  —  Ta'(A))-2  will  be  clear  later.  The  function  Iqv ( A )  is  called  the  gen¬ 
eralized  cross-validation  function  for  the  data  f,  and  the  matrix 

J(\)  =  BkthM(X)BlhW  (91) 

in  Eq.  (90)  is  called  the  influence  matrix.  (See  [3-5,8]). 

Let  A  =  A/v  be  the  A  that  minimizes  the  generalized  cross-validation  function  in  the 
seme  that 


Kv(A.v)  =  nnnVjv(A).  (92) 

In  the  following  we  will  show  that  for  large  data  samples  (i.e.  when  N  is  a  very  large 
value).  A  a-  can  be  used  as  a  good  estimate  of  the  optimal  smoothing  parameter  A*v  in 
the  >ense  that  the  error  in  approximating  the  actual  signal  {/*(r,)}  by  using  the  optimal 

least-squares  estimator  S*(t\  Aa  )  is  asymptotically  the  same  as  that  by  using  S*(t\  A*v) 
with  the  optimal  \*N ,  namely: 


lim 

/V  — oo 


£Tt(\N) 

£Tf(\*N) 


=  1. 


(93) 


(Note  that  the  limit  actually  exists.) 

Before  we  attempt  to  prove  Eq.  (93),  we  need  a  better  understanding  of  £Tf(\). 
First,  by  using  the  notations 


s°(A) 


5*(r,;  A)' 
5*(rN;A) 


and 


'/*(r,)‘ 

f*(TN) 


(94) 


(95) 


we  have,  from  Eq.  (86), 


-VTf(  A)  =  ^,rt(S*(r„  A)  -  /*(r,))2 

1=1 

=  (>°(  A)  -  f)rir;.s°(A)  -  f‘). 

But  since  s°(  A )  =  Dhj,  c*(  A )  =  Bt,t,  M(  A  )Z?^  fi  U  f  =  J  (  A)  f  [sec  F.qs.  (  84  l  a  ml  (  9 1  )  j,  we 
have* 


25 


Chui 


-VTf(A)  =  (J(A)f  -  VfW{J(  Ajf  -  f*) 

=  [j(A)(r  +  e)  -  r]Tu'[j(A)(r  +  s)  -  ri 
=  [J(A)P  -r]rU'[J(A)f*  -  P] 

+  crjT(\)W[j(\)r  -  r]  +  [j{ Ajr  -  r]ru\7(A)c 

+  £TJT(\)]VJ(\)c 

whore  c  =  [t i  .  .  .c»]r.  Hence,  it  follows  from  the  assumptions  £(-,)  =  0  and  £(e,Sj 
<7*b,j  on  the  noise  that 

£Tr(\)  =  l[j( Ajr  -  r]rn'[j(A}r  -  pj 


+  -  Tracc[JT(A)irS2J(A)] 

=  rf.(A)+  f  Trace[^' 

^1=1  ' 


where  E2  =  diag(cr2,. . .  ,  cr2)  has  been  combined  with  the  weight  matrix  W  by  using 


the  choice  of  W  in  Eq.  (49)  to  yield  the  constant  factor  (  al  2)  1  of  the  tract',  and 

i= 1 

Tf  differs  from  Tf  in  Eq.  (86)  in  that  the  optimal  estimator  S*( r,,  A)  in  Eq.  (86)  must 
now  be  replaced  by  the  corresponding  one  that  estimates  the  noiseless  signal  f*.  On  the 
other  hand,  by  the  same  derivation  as  above,  we  also  have 


£Tr(\)  -  £Tf(\)  +  —  ^Trace  TEE2  -  2  Trace  WT?J{ A) 
/  tL 

=  £Tf( A)  +  /  J2  ^  )  f‘v  ~  2  Trace  ^a)]; 

^  i—i  ' 


ETf( A)  =  ^1  -  -^r  Trace  J( A))  £V\(X) 


+  ^/rf2)  [2  Trace  J{\)  —  -V], 


In  the  following  estimates,  in  order  to  clarify  the  presentation,  we  will  use  the  abbrevia¬ 
tions: 


a  =  —  Trace  J{  A  ).  b~  Trace[j7 1  (  A  )J(  A  )’.  r 


Hence,  by  using  Eqs.  (98)  and  (100),  we  liave 
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£Tf( A)  -  £VW(X)  +  cN 
£T{ 

Tf.(A)  +  be  -  (1  -  a)_2[Tf. ( A)  +  be  -  2acN  +  cN ]  +  c.Y 

Tf.(A)  +  fee 

—2a  +  a2  t  (2a cN  —  cN)  +  (cN  —  2acN  +  ca2N) 

(1  -  a)2  +  (Tf-  ( A)  +  bc)(  1  —  a)2 


=  a  —  2  + 


Tf.  +  be  1  (1  —  a)2 


By  setting 


«.vO)  =  -r^Lfp-aH- Av|0 

( 1  —  a  r  o 


‘  v  '  (1  -  a)2  L' 

and  observing  that  Tf.  (A)  >  0.  we  have 

£T{(X)-£Vn(X)  +  cN 
£Tf(\) 


<MA). 


(103) 


(104) 


Now,  by  the  definition  of  A  =  Aw,  we  have  Vjv(A)  <  Vjv(A*),  so  that  £V\(X)  <  £V\(  A* ). 
and  hence  Eq.  ( 104)  implies: 


(£Tf(A))(l  -  6n(A))  <  EFn(A)  -  cN 

<  £Vn(X*)~cN 
<(£Tf(X'))(l  +  6N(X)). 


(105) 


In  addition,  by  the  definition  of  A*  =  X*N,  we  have  £Tf(A*)  <  £" Tf ( A ) .  Therefore,  it 
follows  that 


1  <  1  +  <5;v(A>v) 

£Tf(X*N)  1  —  6n(Xn) 


(106) 


These  inequalities  will  yield  the  limit  result  in  Eq.  (93)  once  we  establish  6 ,\- ( A \  )  — +  0: 
but  in  view  of  its  definition  in  Eq.  (103)  and  the  definition  of  a  in  Eq.  (101).  we  in¬ 
deed  have  6\(X)  — >  0  uniformly  in  A  >  0  if  we  can  show  that  Trace  J{  X)  is  uniformly 
bounded.  This  fact  will  be  verified  in  the  next  section.  We  emphasize  once  again  that 
Eq.  (93)  assures  us  that  A  =  A t\  provides  a  good  estimate  for  the  optimal  smoothing 
parameter  A*  =  A^,. 


ALGORITHMS  FOR  NEAR-OPTIMAL  ESTIMATION 


In  estimation  of  the  signal  onset  of  an  underwater  acoustic  signal  f*  =  {/'*(—,)} 

from  the  noisy  data  f  =  {/, }  where  /,  =  /*(r,)  +  e,,  i  =  1 . -V.  by  applying  the 

method  of  penalized  least-squares  regression  as  discussed  in  the  section  entitled  Penal 
ized  Least -Squares  Estimation,  one  of  the  main  difficulties  is  the  determination  of  the 
optimal  smoothing  parameter  A.  For  a  large  quantity  >Y  of  data  sample-,  a  good  esti 
mate  of  the  optimal  A  =  A*  is  given  by  A  =  A.\  which  minimizes  the  generalized  cros- 
validation  function  l  .v(A)  defined  in  Eq.  (SO).  Since  the  error  functional  7f  i  \ )  can  be 


Chili 


easily  computed  for  each  fixed  value  of  A  (and  of  h)  by  first  computing  the  psuedoinvers< 
M(  A)  in  Eq.  (Si),  the  only  factor  of  I'\(A)  that  remains  to  be  computed  is  the  trace  of 
the  influence  matrix  J{  A)  defined  in  Eq.  (91).  In  the  following,  we  will  simplify  I'v(A) 

so  that  7f(  A)  and  Trace  J(  A)  can  be  computed  simultaneously,  so  that  A  =  Ay  can  be 
determined  more  efficiently. 


Two  Formulations  of  the  Generalized  Cross-Validation  Function 


We  will  first  write  the  generalized  cross-validation  function  l'v(A)  in  a  more  man¬ 
ageable  form.  To  do  so,  we  need  some  notations.  First,  we  can  assume,  without  loss  of 
generality,  that  each  weight  w,  is  non-zero,  since  the  data  information  /,  for  the  zero 
weight  can  be  dropped.  Hence,  we  may  use  the  notations: 


■ 

0 

.  = 

0 

0 

0 

WN  J 

u'\~2  - 

The  positive  definite  banded  symmetric  matrix  Bk  h  defined  in  Eqs.  (71)  and  (GS)  can 
also  be  treated  in  the  same  fashion  by  first  diagonalizing  the  matrix  and  then  taking 
the  positive  square  roots  of  its  eigenvalues.  Hence.  (Bk  ft)]/2  is  also  positive  definite  and 
symmetric  and  ( Bl{  h)l^2(Bk  /()'/2  =  Bk  h.  Recall  that  by  using  Eq.  (84),  the  influence 
matrix  J{\)  can  be  transformed  to  a  nonnegative  definite  symmetric  matrix 


(1081 


where  the  matrix  B{1  h  should  be  replaced  by  the  identity  matrix  if  the  modification  Eq. 

(77)  of  the  normal  equation  in  Eq.  (7G)  is  to  be  studied  (to  simplify  the  computational 
procedure).  In  any  case,  by  setting 


we  may  write 


WiBk.h 


for  Eq.  (7G) 
for  Eq.  (77) 


J(  A)  -  A(  A'7  A  +  A  A/ )+ X  1 . 


(109) 


Mill! 


[We  remark  that  the  matrix  h(k  k  may  be  singular.  In  this  case,  we  recommend  the  model 
in  Eq.  (77)  over  the  model  in  Eq.  ( 7  G ) .  ]  On  the  other  hand,  by  setting 

f=W7f  ,111 


and  noting 


Trace  J[  A  )  —  1 1  ace  J (  A  ) 


we  may  write 
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Kv  ( A )  = 


or  equivalently. 


—  i)TW(J  ( A  j  f  —  f) 
(1  -  jj  Trace  J{  A))2 
T||H.4(J(A)-/)f||^ 

^(Trace[/-  J(\)})*' 


^(Traee[/-  J(\)})2 


(113) 


(114) 


where  the  non-weight ed  (2  norm  [i.e.  j|(.ri, .  . . , .ryv)||22  =  rf  +  •  ■  ■  +  .r2v]  is  used.  (Xote  the 
similarity  with  the  GCY  function  studied  in  Golub,  Heath,  and  Wahba  [S]).  Let  ,Sj  > 
•••'>$(>  S(+i  =  •  •  •  =  s.v  -  0  be  the  singular  values  of  A":  that  is. 


•S1  >  •  •  •  >  1  —  •  •  •  —  -SA’  —  0 


(115) 


are  the  eigenvalues  of  X 1 X.  Then  by  using  the  normalized  eigenvectors  of  XTX  corre¬ 
sponding  to  these  eigenvalues,  we  can  form  two  unitary  matrices  V  and  V  (i.e.  UU1 
U7  U  =  /v  and  VV7  =  V'rV'  =  In+k-i  where  7,v  and  In+k- i  are  identity  matrices)  such 
that 


X  =  UTV‘ 


(116) 


where 


Note  that  .Sf+I  =  ■  •  •  =  „s,v  =  0  and 


r+  = 


1 1 1 5  i 


(See  [1],  pages  1S-22).  H('tice, 
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X  1  X  +  X  XI )+  =  V 


.f+A.v 


O 


o 


»i+\N 


o 


xv  o 


O 

In  addition,  wo  hnvo.  from  Eqs.  (110).  (11G).  and  (119). 

o 

■’1  '  — 

J(X)  =  £' 


O  XV 


O 


d  +  A-V 


o 


:  O 
^  OJ 


v 


•T 


and 


r-J(\)  =  V 


A  V 


*  ?  4-  A V 


O 


o 


A.V 
+  A.V 


o 


Thoroforo.  by  using  f hr  transformation: 

/. 

fx 

dn’ir  f  is  flic  o  'ginal  noisy  data,  wo  havo 


o 


o 


o 


r7f  =  rTir*f 


V 


i=  i 


+  XX 


r-T 


(119) 


(.120) 


(121) 


(122) 


123 ) 
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In  addition,  it  is  clear  that  if  D  is  diagonal  and  U  is  unitary,  then  the  eigenvalues  of 
U DUT  are  the  same  as  those  of  D ,  namely:  the  diagonal  elements  of  D.  Since  the  trace 
of  a  square  matrix  is  the  sum  of  all  its  eigenvalues,  we  have,  from  Eq.  (121).  that 


( 

Trace(/~  J(A))  =  V 
1=1 


A  N 

*1  +  A  N 


+  N-(. 


Putting  Eqs.  (123)  and  (124)  into  Eq.  (114),  we  have 


(124) 


™>-£(^)y(£[£ 

1=1  *  /  \  L 


^  sj  +  NX 

i=i  1 


+ 


n  -e 

~nT 


(125) 


This  formula  is  also  similar  to  Eq.  (2.3)  given  in  Ref.  S. 

Both  formulations  of  Eqs.  (114)  and  (125)  of  I'y(A)  using  the  two  different  trans¬ 
formations  of  the  noisy  data  f  =  and  f  =  UT\V*  f,  respectively,  will  be  useful  for 

computing  the  near  optimal  smoothing  parameter  A  =  A/v. 

For  the  moment,  let  us  apply  Eq.  (120)  to  study  the  truth  of  Eq.  (93).  From  Eqs. 
(106)  and  (108),  we  have  already  noted  that  Eq.  (93)  is  a  consequence  of  the  uniform 
boundedness  of  Trace  J( A)  as  a  function  of  N  for  A  >  Ao  >  0  where  Aq  is  arbitrary.  By 
Eq.  (120),  we  have 


t 

Trace  J (X)  = 

t=i 


s]  +  A N  ' 


(126) 


Hence,  for  any  fixed  A0  >  0,  if  we  can  show  that  s “j  <  M  for  all  i  =  1, . . .  and  all  N, 
then 


0  <  Trace  J( A)  <  ^  <  oo  (127) 

for  all  large  N.  However,  since  are  the  eigenvalues  of  X  r X ,  it  follows  from  Gersch- 
gorin's  theorem,  that 

.v 

max^-  <  max)  ja,,|  (128) 

i  i 

j=i 

where  X1  X  —  [a ,yj.  The  boundedness  of  the  sum  of  the  absolute  values  along  each  row 
can  be  verified  by  using  the  properties  of  5-splines  and  Eq.  (109):  this  is  easily  seen  if 

the  modified  normal  equation  in  Eq.  (77)  is  used  since  X1  X  =  D^WD^t,  is  banded 
and  all  nonzero  entries  are  positive.  A  different  but  much  more  complicated  proof  can  be 
found  in  Refs.  5  and  9. 

Formulation  of  the  Error  Functional 

In  determining  the  time- of- arrival  =  <7  --  nh*.  we  must  determine  h*  which  is 

the  minimum  number  in  the  set  H  =  {/;},  where  each  h  minimizes  the  error  functional 
Ei(li.X)  for  e'ach  fixed  value  of  A  >  0.  Rerall  that  El(Ii.X)  is  the-  error  of  optimal  pe¬ 
nalizes!  estimation  of  any  given  noisy  data  f  from  the  spline  space-  Sq  „  as  define-d  in  Eq. 
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(65).  We  must  now  give  a  computationally  feasible  formulation  of  Eh{h ,  A).  Our  aim  is 

to  arrive  at  a  formula  which  allows  us  to  determine  h  without  getting  into  the  trouble  of 
computing  c*  that  uniquely  determine  the  last  estimator  S £  =  S^(-;  A.  h). 

From  Eq.  (76)  or  Eq.  (77)  and  by  using  Eqs.  (109)  and  (111),  we  have 


c*  = 

=  +  N\I)*Xrt 

where  D ^  h  is  to  be  replaced  by  /  if  Eq.  (77)  is  used.  Hence,  we  have 


c*'  Bt Mc*  =  fT X[(XT X  +  N \I)+Y X 1  f, 
so  that  an  application  of  Eqs.  ( 116)  and  (119)  yields 


+  12  \rT  \ 


(129) 


(130) 


s\  +  A/V 


c*rJ52  lC* 


s*  +  \N 


:  O- 


s2  +  A/V 


(131) 


s2  +  AA' 


4^  Vs?  +  7VA/  ' 

1  =  1 

Again,  Z?j(  is  to  be  replaced  by  /  if  Eq.  (77)  is  used.  Consequently,  by  using  the  for¬ 
mulation  in  Eq.  (162)  where  Sk(t)  and  c  are  to  be  replaced  by  S£(/)  and  c*  respectively, 
and  following  the  same  argument  as  the  one  that  yields  Eq.  (123).  we  have 

Ek(h,  A)  =  Gh(Sl  A) 

=  A|l-Y(-V  '_V  +A.V/)+.Y'f-f||;,  +  Xc-'B'Ij.c'  11321 


4£hrW)^  +  -'Eb 


;  +  'v 


“  V 
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which  can  be  simplified  to  be 


e 

Ek(h,  A)  =  A  y " 

i=i 


s*  +  N\' 


(133) 


An  Algorithm  for  Determining  the  Signal  Onset:  Algorithm  II 


We  are  now  ready  to  give  an  algorithm  for  estimating  the  signal  onset  of  an  under¬ 
water  acoustic  signal  f  =  {/,}  where  f,  =  +  t,  with  £s,  =  0  and  £(t,c;)  =  cr2lbl y 

First,  we  must  choose  the  interval  [0,rf]  for  signal  measurement.  For  more  signal 
information,  we  should  choose  a  large  enough  d  that  by  inspection  the  signal  onset  to  to 
be  determined  should  lie  in  the  first  half  of  the  window;  that  is,  0  <  to  <  d/2. 

Let  jY  be  the  number  of  samples  taken  at  rj , . . . ,  r/v  where  0  <  rj  <  ■  ■  •  <  t  y  <  d. 
We  could  use  equally  spaced  r,’s  say: 

Ti  =  Jf,  i  =  l,...,N.  (134) 

Usually.  N  is  a  fairly  large  number,  say  iV  =  50  to  100.  In  addition  to  the  sampling 
quantity  N ,  we  must  also  decide  on  the  number  n  of  knots  of  the  spline  space,  the  dis¬ 
cretization  of  the  h  values  in  determining  H  =  {/i},  and  the  initial  choice  as  well  as  rate 
of  decrease  in  estimating  A  from  Vm(\). 

We  will  choose  10  <  n  <  50,  h  =  hJ ,  where 

h}  =  -  -  7777-1  j  =  0, . . . ,  M,  (135) 

n  2  M  n 

( M  being  the  number  of  discretization  samples  for  finding  h  and  typically  we  may 
choose  M  cA),  where  h°  and  h correspond  to  the  signal  onset  at  0  and  d/2,  respec¬ 
tively.  For  a  global  search  of  A  to  minimize  V/v(  A),  we  pick  an  upper  bound  A0  and  work 
our  way  down  to  arrive  at  an  optimal  A.  Typically  this  upper  bound  A0  should  depend 
on  the  signal-to-noise  ratio  “5/A”;  the  smaller  this  ratio  the  larger  A0  should  be.  since 
the  signal  {/,}  would  be  fairly  noisy.  In  the  engineering  literature,  ilS/N"  may  be  de¬ 
fined  as  the  ratio  of  the  expectation  of  the  signal  over  the  standard  derivation  a,  of  the 
noise.  For  computational  simplicity,  we  may  define 


“S/Ar"  = 


;V 

yI>,2 

i=i 


a£(£/.)2 

i=i _ 

N 

jjZ-l 

1=1 


(136) 


In  underwater  acoustic,  this  number,  which  is  also  called  the  signal-power/noise- power 
ratio,  is  usually  in  the  neighborhood  of  40  db.  We  recommend  choosing 

A0  =  c/“S/.V”  i  137  i 

since  the  noisier  the  signal,  the  larger  the  smoothing  parameter  U  required.  I  In'  < •  i i *  > : < •  <  ■ 
of  the  constant  r  must  depend  on  the  experiment.  Following  Craven  and  Wahha  '  we 
pick 
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AJ  =  AJ_1 10~^ .  j  =  1.2 . 

for  the  search  of  A  =  A,v-  (Note  that  \J  — ►  0  as  j  — »  oc,  and  A"  is  suppos('d  to  be  an 
upper  bound  of  A.) 

The  computational  procedure  can  be  summarized  as  follows: 


Algorithm  II  (Spline  Estimation  of  Signal  Onset  in  Noisy  Environment) 


(1°) 


(2°) 

(3°) 

(4°) 


(0°) 


(S°) 

(9°) 

(10°) 

(11°) 


f  12°) 
I  13°) 


Det  ermine  of,  i  -  1  and  apply  the  formula  in  Eq.  (49)  to  find  the  weights 

ir,.  i  =  1, ....  Ad  Form  the  weight  matrix 


and  its  square  root 


o 


H'X 


( 138) 


’  x/iET 

O  ' 

.  o 

\furN  _ 

(139) 


Input  values  h  =  hJ .  j  =  0 . M.  [See  Eq.  (135)]. 

For  each  h  —  h] ,  compute  in  Eq.  (GS)  to  form  the  positive  definite  symmetric 

matrix  Df.  h  in  Eq.  (71).  For  k  =  1.3.  this  has  been  done.  [See  Eq.  (SO)  for  k  =  3j. 
Compute  the  inverst'  of  the  positive  square-root  ( D ^  /()-1/*  of  Df  ,  . 


[Note  that  if  the  modification  Eq.  (77)  of  Eq.  (7G)  is  used,  skip  (3°)  and  (4°)  and 
proceed  to  (5°).  replacing  ( D jj?  h)~ i  by  /.] 


Compute  Dk,u  in  Eq.  (74). 

Compute  X  =  Xh  in  Etp  (109). 

Determine  the  SVD  (singular  value  decomposition)  of  Ad  I  T]’1  where  V  —  L),  and 
V'  =  l>  are  unitary  matrices  and  T  =  T/,  dej>ends  on  h.  [Set'  Eqs.  (11G)  and  (  117)]. 

Compute  the  transformed  data  f  =  [’Mid  f. 

Input  values  A  =  X} , \J  =  A-7  —  1 10  ,  j  -  1.2....,  by  using  an  appropriate  e  in  E< } . 

(137)  to  givt'  an  initial  value  A°. 

Compute  I'v(A)  in  Eq.  (GG)  by  using  (6°)  and  (7°)  for  A  =  AJ  in  (9°). 

Determine  A  =  A,v  among  {A7}  so  that 


Vs(\)=r  min  <139- 

7  =  0. 1  . . . 

Compute  Eh(h.  A)  in  Etp  (133)  for  h  —  h1 .  j  =  0 . M. 

Determine  h  among  hJ  so  that 

E|(//.A)=  min  Er(/d,  A).  (Mil 

j  =  0.  ...At 
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(14°)  Select  the  smallest  h*  from  H  —  {/i}  in  (13°). 
(15°)  Evaluate  the  signal  onset 


tQ  =  d  —  7  ih* . 


(141) 


(This  is  the  required  answer). 

To  plot  the  spline  curve,  proceed  to  the  next  two  steps. 
(16°)  Compute  c*  by  using  the  formula 


c*  = 


sJ  +  A  N 


o 


and  write  c*  —  [clt+1 , . . . ,  c*_j] 
(17°)  Plot  the  spline  curve  of 


n  —  1 


O 


O 


n^  +  XN 


Sk(t)=  Y 

>=-*+1 


f 


o 


O-l 


(142) 


(143) 


which  gives  a  smooth  approxiinant  of  the  noisy  signal  with  the  "optimal”  signal  on¬ 
set  t0. 


Chui 


A  flow-chart  of  this  algorithm  is  given  below. 


/  INPUT 
e  >  0,  d,  n,  M,  N,  I, 
arrays  T,  f,  XV,  B°k  h 

e  <—  0  ,L*—( 

3,  K  —  0 

(Note  A) 


h  «—  d/n 
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Note  A: 

e  -  a  given  small  positive  number, 
d  -  an  upper  bound, 

I  -  search  width,  can  be  selected  to  be  (as  suggested  in  the  algorithm), 
n  -  number  of  knots,  can  be  selected  as  large  as  50, 

M  -  search  density,  can  be  selected  as  M  =  cN  for  some  suitable  constant  c. 
N  -  the  number  of  observation  points, 

T  =  [tj  , , . , ,  tn]t  -  ordinate  vector  of  observation  points, 
f  =  [/i , . . . ,  fs)7  -  noisy  data  set  at  sample  T, 


-tCi 

°1 

■o 

U’N 

the  matrix  of  weights  with  square  root  W1/2 , 


Z?1/2  <—  (B%  h)  inverse  of  the  positive  square-root  of  Z?£  h 


Note  B:  For  k  =  3  (cubic  splines),  B  =  ZZ*  /,  can  be  easily  computed  as  follows: 

Set  B  = 

( 1°)  For  i  =  1, . . . ,  TV;  j  =  3, . . . ,  n  +  2,  do  the  following: 

Let  ZZ  —  n  -j-  3  —  j  —  j^(d  —  t,). 

If  0  <  H  <  1  then  BtJ  =  |ZZ3, 

if  1  <  H  <  2  then  Bij  |(2  -  ZZ)3  +  (2  —  ZZ)2(ZZ  -  1)  :  2(2  -  ZZ)(ZZ  -  l)2  +  §  (ZZ  -  1)\ 

if  2  <  H  <  3  then  B.y  =  |(3  -  ZZ)3  +  2(3  -  ZZ)2(ZZ  -  2)  +  (3  -  ZZ)(ZZ  -  2)2  +  ±(ZZ  -  2)3. 

if  3  <  H  <  4  then  ZZ.y  =  |(4  -  ZZ)3; 

otherwise  BtJ  =0. 

(2°)  For  i  =  1 , ,N;j  =  2,  do  the  following: 

Let  H  =  \(tt  -  t0). 

If  0  <  H  <  1  then  Bl}  =  |(1  -  ZZ)ZZ2  +  ^ZZ3, 

if  1  <  ZZ  <  2  then  ZZ,,  =  £(2  -  ZZ)3  +  2(2  -  Z/)2(ZZ  -  1)  +  (2  -  ZZ)(ZZ  -  l)2  +  \(H  -  l)3 

if  2  <  ZZ  <  3  then  ZZU  =  |(3  -  ZZ)3; 

otherwise  BtJ  =  0. 

(3°)  For  i  =  1, . . . ,  A,  j  =  1,  do  the  following: 

Let  ZZ  =  £(t,  -  t0)- 

If  0  <  ZZ  <  1,  then  ZZ,,  =  3(1  -  ZZ)2ZZ  +  f  (1  -  ZZ)ZZ2  +  |ZZ3, 
if  1  <  ZZ  <2  then  ZZ0  =  J(2  -  ZZ)3: 
otherwise,  ZZ,,  =  0. 

Note  C:  Determine  SVD  of  A'  =  ZTV'7;  where  Zr  and  V  are  unitary  matiicr.-  and  F  is 
given  in  Eq.  (117). 

Note  D:  The  matrix  formulations  o/V',v(A)  and  Ei,(h.  A): 

Write 
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and 


Then 


and 


iVfTffrf 

V)v(A)  = - ^3 - 

;  (eT(rr7’)I/2e)2 

Ek(h,  A)  =  Afr(frT)1/2f. 


(144) 


(145) 


(146) 

(147) 


Another  Algorithm  for  Determining  the  Signal  Onset:  Algorithm  III 

Instead  of  finding  the  SVD  of  the  N  x(n  +  k—l)  matrix  X  as  in  Eq.  (116)  to  obtain 
a  diagonization  of  XTX  +  X  XI,  it  is  more  efficient  to  tridiagonalize  the  symmetric  square 
matrix  XX  r  (see  Golub  and  Van  Loan  [10]), 

XX  t  =  UTUt  ( 14S) 

where  U  is  unitary  and  T  is  nonnegative  symmetric  and  tridiagonal.  Note  that  U  can  be 
stored  in  factored  form  in  the  strict  lower  triangle  of  XX1  (see  [10].  pages  276-277,  and 
Gu  and  Wahba  [9],  page  8).  In  addition,  since  T  is  tridiagonal,  it  is  easy  to  find  the  CD 
(Cholesky  decomposition)  of  N XI  +  T,  namely: 

ArA/  +  f  =  rTr  (149) 

where 

‘«i  b] 

Y=  ''  '■  O  (150) 

O  «/V-l  ^.v-l 

«.V  - 

Note  that  since  T  is  symmetric  and  nonnegative  definite,  XX I  +  T  i>  invertible  for  all 
A  >  0.  so  that.  (i\. . . .  ,  (i ,\  are  all  nonzero. 
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We  can  now  study  the  generalized  cross-validation  function  V)v(A)  as  expressed  in 
the  form  Eq.  (114)  as  follows.  From  Eqs.  (121)  and  (116).  consecutively,  we  have 

/  _  J( A)  =  XNU(TTt  +  \NI)+Ut 

=  \NU(UtXXtU  +  A  NI)+Ut  ( 15 

=  \N(XXt  +  \NI)+ . 

so  that  Eq.  (114)  becomes 

„  ||(.V.Y3'  +  A.V/)+f||^ 


(152) 


We  now  use  the  unitary  matrix  U  in  Eq.  (148)  to  make  another  data  transformation: 


x  =  UTf  =  UTW${,  (153) 

so  that  by  using  Eqs.  (148)  and  (149),  Eq.  (152)  can  be  written  as 

,,  \\(YTY)^rtl 

VnW~J[ WW'  (154) 

As  mentioned  above,  for  A  >  0  (which  will  be  assumed  in  the  following  discussion),  F7  F 
is  invertible,  so  that  (FrF)+  =  (FrF)-1,  and  F/v(A)  takes  on  the  form: 


(154) 


v  (\\  -  l^xr(FrF)-2x 

NK  )~  (jr  Trace(F7’F)“1  )2 


(155) 


To  discuss  the  computational  procedure  in  minimizing  Eq.  (155),  let  us  use  the  no¬ 
tation 

(Ft)_1  =  ...yj 

where  y^_i+1  denotes  the  ith  column  vector  of  (F-1)r.  This  yields 


I  =  (Y'1)tYt  =  [y*  ...yj 


(157) 


which  is  equivalent  to 


bN  aN 


yi  =  —  e,v 

U  V 


y2  = - (ew_,  -  y, 

a\~i 


(15S) 


y.v  =  —  (e,  -  ftiy.v_,  )• 
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where 


with  the  value  1  in  the  ilh  component,  denotes  the  ith  standard  unit  vector  in  1R  N  . 

Next,  since  ( V  ~ 1  )T  is  a  lower  triangular  matrix,  its  (i  +  l)st  column  y,v_,  is  orthogonal 
to  e,.  Hence,  Eq.  (158)  yields: 


|y.ll#>  =  aN 


lly^ltf2  -  aN-i(i  +  &N-illy. lit2' 


(ICO) 


l  IlyJlp  =  «t  2(i  +  ^llyyv-,11?2)- 

The  importance  of  Eq.  (160)  is  that  it  yields  an  efficient  computational  procedure  of  the 
denomination  of  in  Eq.  (155)  since 


1  Tracer)-1  | 


Computation  of  the  numerator  of  V)v(A)  in  Eq.  (155)  simply  involves  four  matrix-vector 
multiplications  followed  by  a  vector-vector  multiplication: 

f1  :=  W?  f 


=  ulv 


f3:=(FT)-if2  (1G2) 

f4  :=  y - 1  f 3 

ixr(v_7'y)-,x  =  if*  ■  f*  =  i(f*)Tf*. 

Here,  VV’’*  is  a  diagonal  matrix,  and  (Y'T)~l  and  V-1  are  triangular  matrices  given  bv 
Eq.  (158). 

The  final  formula  we  need  is  one  concerning  the  error  functional  £’*(5,  A).  For  com¬ 
puting  this  error  functional  in  terms  of  the  information  given  by  the  new  decomposition 

(i.e.  V,  7\  etc.),  we  also  need  the  pseudo- inverse  V0+  of  V)  when* 

f  =  Yory0.  163) 

[Here,  it  should  be  noted  that  F’o  can  be  identified  as  }’  in  Eq.  (150)  by  taking  the  limit 
as  A  — +  0.  See  Eqs.  (149)  and  (150)].  So.  Vo  may  be  singular.  Fortunately,  the  hi 
diagonal  matrix  V0  is  independent  of  A.  so  that  for  each  value  of  h.h  —  hJ .  the  p>*  udo- 
inverse  can  be  computed  and  stored  without  any  reference  to  the  changes  of  A.  F><r  nota- 
tional  consistency,  we  write 
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[y°  y° 

l*7  /V  J  TV  -f  1 


(164) 


where  y^  is  the  ?,h  column  vector  of  (V+)T- 

We  are  now  ready  to  derive  a  formula  for  Ek ( /> ,  A ) .  For  this  purpose,  we  need  tin' 
results  from  the  SVD  of  A  in  Eq.  ( 1 1G )  as  intermediate  steps.  (Of  course,  the  final  for¬ 
mula  must  not  depend  on  the  SVD  of  A'.)  First,  observe  that 


[/  -  N\(XXr  +  7VA/)-1](A'Ar)+[7  -  XX(XXT  +  ArA)_I] 

=  [i-  nx(XXt  +  N\iyl}U{rrT)+uT[i  -  xx(xxr  +  aa)_1] 

=  F(rr7’  +  Nxir'rrTuru(rrT)+uTU(rrT  +  xxir'rr't'u1'  (i05) 
=  u(  rrT  +  NXi)~2rrTuT  =  ur(rrr  +  xxi)~2rTuT 

=  X(XTX  +  NXI)~2Xt. 


Hence,  it  follows  that 


c*1  B°k  hc*  =  tTX{XTX  +  NXI)~2XTf 

=  fr[7  -  NX(XXt  +  7VA7r‘](AAr)+[7  -  XX(XXr  +  AA)~!jf 
=  fr[7-  NX(UTUt  +  NXI)~1](UTUt)+[I  -  NX(UTUt  +  NXI)~'}f  (166) 
=  fTU[I  -  NX(T  +  NXI)~l)f+[I  -  NX(T  +  iVA7)_1]f7Tf 
=  xr[7  -  iVA(rry)-1](ir0Tr0)+[7  -  nx(yty)-1]x, 

where  Eqs.  (148),  (149),  and  (163)  have  been  used.  Hence,  by  using  the  numerator  of 
Eq.  (155)  [see  Eqs.  (114)  and  (152)],  we  have 

&(M>  =  A|l(r7'VT,x||J,  +  A||(v;,+)V-  iVA(y7Vr')x||^.  (i67> 

In  computation,  we  may  write 

Ek{h,l)  =  Xtffi*  +  xHV+flt1  -  a\ti)\THY+)T(t1  -  N\t*)\  (108) 

where  f2  and  f4  have  been  computed  in  Eq.  (162)  and  V0+  is  taken  from  Eq.  (164). 

After  A  =  A  and  h  =  h*  are  determined,  we  may  wish  to  plot  the  spline  curve. 
Hence,  a  formula  for  the  77-spline  coefficients  c*  =  [elt+1  .  .  .  c*_,]r  is  required.  Since 

A  and  h  are  already  fixed  (to  be  A  and  h *  respectively),  computational  efficiency  is  not 
as  important.  In  any  case,  we  will  use  as  much  known  information  as  possible.  However, 
the  SVD  of  A*  for  determining  A+  does  not  seem  to  be  avoidable.  First  note  that  since 
X  =  UTV1 ,  we  have 


A + A  (A 7  A  +  AA7)-'A7 

=  vr+r(rrr  +  AAfr'rr7 
=  v(rTr  +  xxir'Tu1 
=  (XTx  +  xxiy'x1. 


( 169) 
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It  should  be  emphasized  that  in  the  above  derivation  r+T  is  not  necessarily  the  identity 
matrix.  Fortunately,  the  zero  block  in  its  lower  right  corner  matches  that  of  T  which  also 
occurs  in  the  same  expression  above.  Next,  as  before,  we  also  have 


X(XrX  +  XX  I)~lXrf 
=  [I-  XX(XXr  +  XXI)-l}{. 


( 170) 


Hence.  from  E<p  (7G)  or  Eq.  (77)  and  by  using  Eqs.  (109)  and  (111),  we  have 

c*  =  (Dh,  rUxTx  +  xxir1  xTi 

=  (B°kj,rl*X+ X(XrX  +  X  XI)~l  X 1  { 

=  (Bl,,rxx+[I  -  XX(XXT  +  XXI)~]}{  (171) 

=  (B°k  h)-XX+U[I  -  XX (T  +  NXI)~l]Urf 

=  (Blh)~XX+U[I  -  XX(YrY)-'}x. 

In  computation,  by  using  Eq.  (162),  this  is  equivalent  to 

c*  =  ( B°k  h  )  ~  XX  +  U  { f'2  -  iVAf4  } .  (172) 

Of  course,  when  the  modified  penalized  factor  Xc1  c  in  place  of  Ac7  (B°k,  h)c  is  used, 
then  we  have 


c*  =  X+U{(2  -  iVAf1}.  (173) 

Summarizing  the  above  discussions,  we  have  the  following  computational  procedure: 

Algorithm  III  ( Spline  Estimation  of  Signal  Onset  in  Noisy  Environment  II) 

(1°)  Same  as  (1°)  in  Algorithm  II. 

(2°)  Same  as  (2°)  in  Algorithm  II. 

(3°)  Same  as  (3°)  in  Algorithm  II. 

(4°)  Same  as  (4°)  in  Algorithm  II. 

[Note  that  if  the  modification  in  Eq.  (77)  of  Eq.  (76)  is  used,  skip  (3°)  and  (4°)  and 
proceed  to  (5°),  replacing  ( B £  h)~^  bv  I.] 

(5°)  Same  as  (5°)  in  Algorithm  II. 

(6°)  Same  as  (6°)  in  Algorithm  II. 

(7°)  Compute'  tridiagonal  matrix  T  and  unitary  matrix  V  such  that  XX1  =  VTV  1 . 

(8°)  Same  as  (9°)  in  Algorithm  II. 

(9°)  Compute  CD  (Choleskv  decomposition): 

XXI  +  T  =  I  1  V,  V  as  given  in  Eq.  (150).  i  174) 

(10°)  Compute  ijx . y,v  by  using  Eq.  (158). 

(11°)  Compute  Trace(V7  Y)~l  in  Eq.  (161)  by  using  Eq.  (160)  and  (10°). 

(12°)  Compute  f 1 ,  f 2 .  f 5  and  f4.  and 
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(13°)  Compute  V'v(A)  =  - )  f - 7  bv  using  (11°)  and  (12°). 

( Tr ace(  V  r  Y' )  ~  1  ) 

(14°)  Same  as  (11°)  in  Algorithm  II. 

(15°)  Compute  CD  (Cholesky  deeomposition)  T  =  l’"7}'0  [i.t>.  (9°)  with  A  =  0]  and 
compute  SVD  to  form  y()+  as  in  Eq.  (1G4). 

(1G°)  Compute  £*(/>,  A)  in  Eq.  (1GS)  by  using  A  from  (14°)  for  h  —  hJ.  j  =  0 . \I. 

(17°)  Same  as  (13°)  in  Algorithm  II. 

(IS0)  Same  as  (14°)  in  Algorithm  II. 

(19°)  Same  as  (15°)  in  Algorithm  II. 

(20°)  Compute  A’+  using  A  from  (14°)  and  h*  from  (19°). 

(21°)  Compute  c*  by  using  Eq.  (172)  or  Eq.  (173)  depending  on  the  mathematical  model 
for  penalized  estimation. 

(22°)  Same  as  step  (17°)  in  Algorithm  II 

A  flow-chart  for  this  algorithm  is  given  below. 
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